Algoritmo para matrices tridiagonales

El algoritmo para matrices tridiagonales o algoritmo de Thomas (por Llewellyn Thomas) es un algoritmo del álgebra lineal numérica para resolver matrices tridiagonales de forma eficiente.

Una matriz tridiagonal se corresponde a un sistema de ecuaciones de la forma

$a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}=d_{i},\,\!$

donde $a_{1}=0\,$ y $c_{n}=0\,$. lo que se puede representar matricialmente como

${\begin{bmatrix}{b_{1}}&{c_{1}}&{}&{}&{0}\\{a_{2}}&{b_{2}}&{c_{2}}&{}&{}\\{}&{a_{3}}&{b_{3}}&\ddots &{}\\{}&{}&\ddots &\ddots &{c_{n-1}}\\{0}&{}&{}&{a_{n}}&{b_{n}}\\\end{bmatrix}}{\begin{bmatrix}{x_{1}}\\{x_{2}}\\{x_{3}}\\\vdots \\{x_{n}}\\\end{bmatrix}}={\begin{bmatrix}{d_{1}}\\{d_{2}}\\{d_{3}}\\\vdots \\{d_{n}}\\\end{bmatrix}}.$

Para este tipo de sistemas se puede obtener con este algoritmo una solución con solo $O(n)$ operaciones en vez de las $O(n^{3})$ que requiere la eliminación gaussiana. El algoritmo primero elimina las $a_{i}$ y luego usa una sustitución para obtener la solución.

Este tipo de matrices suelen salir al plantear discretizaciones por métodos de diferencias finitas, volúmenes finitos o elementos finitos de problemas unidimensionales. Algunos de los problemas físicos que se plantean así son la Ecuación de Poisson, la ecuación de calor, la ecuación de onda o la interpolación por splines.

Método

El primer paso del método es modificar los coeficientes como sigue:

$c'_{i}={\begin{cases}{\begin{array}{lcl}{\cfrac {c_{i}}{b_{i}}}&;&i=1\\{\cfrac {c_{i}}{b_{i}-c'_{i-1}a_{i}}}&;&i=2,3,\dots ,n-1\\\end{array}}\end{cases}}\,$

donde se marcan con superíndice ' los nuevos coeficientes.

De igual manera se opera:

$d'_{i}={\begin{cases}{\begin{array}{lcl}{\cfrac {d_{i}}{b_{i}}}&;&i=1\\{\cfrac {d_{i}-d'_{i-1}a_{i}}{b_{i}-c'_{i-1}a_{i}}}&;&i=2,3,\dots ,n.\\\end{array}}\end{cases}}\,$

a lo que se llama barrido hacia adelante. A continuación se obtiene la solución por sustitución hacia atrás:

$x_{n}=d'_{n}\,$

$x_{i}=d'_{i}-c'_{i}x_{i+1}\qquad ;\ i=n-1,n-2,\ldots ,1.$

Implementaciones

Implementación en C

La siguiente función en C resolverá el sistema (aunque sobreescribirá el vector de entrada c en el proceso). Se ha de notar que aquí los subíndices están basados en el cero, es decir $n=0,1,\dots ,N-1$ donde $N$ es el número de ecuaciones:

void solve_tridiagonal_in_place_destructive(float x[], const size_t N, const float a[], const float b[], float c[]) {
    int n;
    
    /**
     *  solves Ax = v where A is a tridiagonal matrix consisting of vectors a, b, c
     *  note that contents of input vector c will be modified, making this a one-time-use function
     *  x[] - initially contains the input vector v, and returns the solution x. indexed from [0, ..., N - 1]
     *  N   - number of equations
     *  a[] - subdiagonal (means it is the diagonal below the main diagonal) -- indexed from [1, ..., N - 1]
     *  b[] - the main diagonal, indexed from [0, ..., N - 1]
     *  c[] - superdiagonal (means it is the diagonal above the main diagonal) -- indexed from [0, ..., N - 2]
     */
    
    c[0] = c[0] / b[0];
    x[0] = x[0] / b[0];
    
    /* loop from 1 to N - 1 inclusive */
    for (n = 1; n < N; n++) {
        float m = 1.0f / (b[n] - a[n] * c[n - 1]);
        if(n < (N-1)) c[n] = c[n] * m;
        x[n] = (x[n] - a[n] * x[n - 1]) * m;
    }
    
    /* loop from N - 2 to 0 inclusive */
    for (n = N - 2; n >= 0; --n) {
        x[n] = x[n] - c[n] * x[n + 1];
    }
}

La siguiente variante preserva el sistema de ecuaciones para reutilizarlo en otras funciones. Se hacen llamadas a la biblioteca para reservar más especiao. Otras variantes usan un puntero a memoria disponible.

void solve_tridiagonal_in_place_reusable(float x[], const size_t N, const float a[], const float b[], const float c[]) {
    size_t n;
    
    /* allocate scratch space */
    float * const cprime = malloc(sizeof(float) * N);
     
    cprime[0] = c[0] / b[0];
    x[0] = x[0] / b[0];
    
    /* loop from 1 to N - 1 inclusive */
    for (n = 1; n < N; n++) {
        float m = 1.0f / (b[n] - a[n] * cprime[n - 1]);
        if(n < (N-1)) cprime[n] = c[n] * m;
        x[n] = (x[n] - a[n] * x[n - 1]) * m;
    }
    
    /* loop from N - 2 to 0 inclusive */
    for (n = N - 2; n > 0; --n)
        x[n] = x[n] - cprime[n] * x[n + 1];
      
    /* free scratch space */
    free(cprime);
}

Implementación en Python

La siguiente implementación usa el lenguaje de programación Python.De nuevo, los subíndices son basados en el cero ($i=0,1,\dots ,n-1$ donde $n$ es el número de incógnitas).

def TDMASolve(a, b, c, d):
    n = len(d)  # número de filas

    # Modifica los coeficientes de la primera fila
    c[0] /= b[0]  # Posible división por cero
    d[0] /= b[0]

    for i in range(1, n):
        ptemp = b[i] - (a[i] * c[i-1])
        c[i] /= ptemp
        d[i] = (d[i] - a[i] * d[i-1])/ptemp

    # Sustitución hacia atrás
    x = [0 for i in range(n)]
    x[-1] = d[-1]

    for i in range(-2, -n-1, -1):
        x[i] = d[i] - c[i] * x[i+1]

    return x

Implementación en Matlab

En Matlab/Octave el algoritmo queda como sigue. Esta vez los vectores están basados en el uno, por lo que $i=1,2,\dots ,n$ con $n$ que el número de incógnitas.

function x = TDMAsolver(a,b,c,d)
%a, b, c are the column vectors for the compressed tridiagonal matrix, d is the right vector
% N is the number of rows
N = length(d);
 
% Modify the first-row coefficients
c(1) = c(1) / b(1); % Division by zero risk.
d(1) = d(1) / b(1); 
 
for n = 2:1:N
    temp = b(n) - a(n) * c(n - 1);
    if (n<N)
        c(n) = c(n) / temp;
    end
    d(n) = (d(n) - a(n) * d(n - 1)) / temp;
end
 
% Now back substitute.
x(N) = d(N);
for n = (N - 1):-1:1
    x(n) = d(n) - c(n) * x(n + 1);
end
end

Implementación en Fortran 90

Fortran usa también una nomenclatura basada en el uno, es decir $i=1,2,\dots ,n$ siendo $n$ el número de incógnitas.

Algunas veces no es deseable que el programa sobreescriba los coeficientes (por ejemplo para resolver diversos sistemas que solo difieren en el término independiente), así que esta implementación mantiene dichos coeficientes.

      subroutine solve_tridiag(a,b,c,d,x,n)
      implicit none
!	a - sub-diagonal (means it is the diagonal below the main diagonal)
!	b - the main diagonal
!	c - sup-diagonal (means it is the diagonal above the main diagonal)
!	d - right part
!	x - the answer
!	n - number of equations

        integer,intent(in) :: n
        real(8),dimension(n),intent(in) :: a,b,c,d
        real(8),dimension(n),intent(out) :: x
        real(8),dimension(n) :: cp,dp
        real(8) :: m
        integer i

! initialize c-prime and d-prime
        cp(1) = c(1)/b(1)
        dp(1) = d(1)/b(1)
! solve for vectors c-prime and d-prime
         do i = 2,n
           m = b(i)-cp(i-1)*a(i)
           cp(i) = c(i)/m
           dp(i) = (d(i)-dp(i-1)*a(i))/m
         enddo
! initialize x
         x(n) = dp(n)
! solve for x from the vectors c-prime and d-prime
        do i = n-1, 1, -1
          x(i) = dp(i)-cp(i)*x(i+1)
        end do

    end subroutine solve_tridiag

Derivación

Algoritmo

Se puede obtener dicho algoritmo usando la eliminación gaussiana de forma genérica. Suponiendo como incógnitas $x_{1},\ldots ,x_{n}$ y con ecuaciones a resolver:

${\begin{aligned}b_{1}x_{1}+c_{1}x_{2}&=d_{1};&i&=1\\a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}&=d_{i};&i&=2,\ldots ,n-1\\a_{n}x_{n-1}+b_{n}x_{n}&=d_{n};&i&=n.\end{aligned}}$

Modificando la segunda ecuación a partir de la primera obtenemos:

$({\mbox{ecuación 2}})\cdot b_{1}-({\mbox{ecuación 1}})\cdot a_{2}$

lo que da:

$(a_{2}x_{1}+b_{2}x_{2}+c_{2}x_{3})b_{1}-(b_{1}x_{1}+c_{1}x_{2})a_{2}=d_{2}b_{1}-d_{1}a_{2}\,$

$(b_{2}b_{1}-c_{1}a_{2})x_{2}+c_{2}b_{1}x_{3}=d_{2}b_{1}-d_{1}a_{2}\,$

y se ha eliminado $x_{1}$ de la segunda ecuación. Si repetimos usando la segunda ecuación en la tercera obtenemos:

$(a_{3}x_{2}+b_{3}x_{3}+c_{3}x_{4})(b_{2}b_{1}-c_{1}a_{2})-((b_{2}b_{1}-c_{1}a_{2})x_{2}+c_{2}b_{1}x_{3})a_{3}=d_{3}(b_{2}b_{1}-c_{1}a_{2})-(d_{2}b_{1}-d_{1}a_{2})a_{3}\,$

$(b_{3}(b_{2}b_{1}-c_{1}a_{2})-c_{2}b_{1}a_{3})x_{3}+c_{3}(b_{2}b_{1}-c_{1}a_{2})x_{4}=d_{3}(b_{2}b_{1}-c_{1}a_{2})-(d_{2}b_{1}-d_{1}a_{2})a_{3}.\,$

Esta vez se ha eliminado $x_{2}$. El procedimiento se puede repetir hasta la fila enésima, dejando cada ecuación con solo dos incógnitas menos la última que solo tiene una incógnita. Entonces es inmediata la resolución de esta y desde ahí la de las anteriores.

Fórmulas para coeficientes

La determinación de los coeficientes en las ecuaciones generales es más complicada. Examinando el procedimiento, se pueden definir de forma recursiva:

${\tilde {a}}_{i}=0\,$

${\tilde {b}}_{1}=b_{1}\,$

${\tilde {b}}_{i}=b_{i}{\tilde {b}}_{i-1}-{\tilde {c}}_{i-1}a_{i}\,$

${\tilde {c}}_{1}=c_{1}\,$

${\tilde {c}}_{i}=c_{i}{\tilde {b}}_{i-1}\,$

${\tilde {d}}_{1}=d_{1}\,$

${\tilde {d}}_{i}=d_{i}{\tilde {b}}_{i-1}-{\tilde {d}}_{i-1}a_{i}.\,$

Para acelerar la resolución ${\tilde {b}}_{i}$ puede ser dividido (si no hay problemas por división por cero) y los nuevos coeficientes, indicados con primas, serán:

$a'_{i}=0\,$

$b'_{i}=1\,$

$c'_{1}={\frac {c_{1}}{b_{1}}}\,$

$c'_{i}={\frac {c_{i}}{b_{i}-c'_{i-1}a_{i}}}\,$

$d'_{1}={\frac {d_{1}}{b_{1}}}\,$

$d'_{i}={\frac {d_{i}-d'_{i-1}a_{i}}{b_{i}-c'_{i-1}a_{i}}}.\,$

Lo que da el siguiente sistema con las mismas incógnitas y los coeficientes definidos en función de los originales:

${\begin{array}{lcl}x_{i}+c'_{i}x_{i+1}=d'_{i}\qquad &;&\ i=1,\ldots ,n-1\\x_{n}=d'_{n}\qquad &;&\ i=n.\\\end{array}}\,$

Donde la última ecuación solo afecta a una incógnita. Resolviéndola podemos resolver la anterior y así sucesivamente:

$x_{n}=d'_{n}\,$

$x_{i}=d'_{i}-c'_{i}x_{i+1}\qquad ;\ i=n-1,n-2,\ldots ,1.$

Variantes

En algunas situaciones, particularmente con condiciones de contorno periódicas, se busca resolver una forma perturbada de la matriz tridiagonal:

${\begin{aligned}a_{1}x_{n}+b_{1}x_{1}+c_{1}x_{2}&=d_{1},\\a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}&=d_{i},\quad \quad i=2,\ldots ,n-1\\a_{n}x_{n-1}+b_{n}x_{n}+c_{n}x_{1}&=d_{n}.\end{aligned}}$

Para este caso se usa la Fórmula Sherman-Morrison para evitar las operaciones adicionales de una eliminación gaussiana y poder seguir usando el algoritmo de Thomas. Este método requiere resolver una versión no cíclica del sistema para la entrada y un vector de corrección y combinar los resultados. Todo ello se puede hacer eficientemente de una vez dado que ambos problemas compartan el barrido hacia delante de la matriz triangular..

En otras situaciones el sistema de ecuaciones puede ser tridiagonal por bloques, con submatrices como elementos formando tres diagonales. Este tipo de problemas suelen darse cuando se quiere extender los métodos de discretización a dos dimensiones y puede verse también como una matriz pentadiagonal. Variantes de este algoritmo se usan para estas situaciones.

El libro Numerical Mathematics de Quarteroni, Sacco y Saleri, recoge una versión modificada para multiplicaciones en vez de divisiones, lo que es útil en ciertas arquitecturas.

Referencias

• Conte, S.D., and deBoor, C. (1972). Elementary Numerical Analysis. McGraw-Hill, New York. ISBN 0070124469. 
• Este artículo incluye texto del artículo Tridiagonal matrix algorithm - TDMA (Thomas algorithm) publicado con licencia GNU en CFD online wiki
• Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). «Section 2.4». Numerical Recipes: The Art of Scientific Computing (3rd edición). New York: Cambridge University Press. ISBN 978-0-521-88068-8. Archivado desde el original el 4 de marzo de 2016. Consultado el 1 de febrero de 2013. 

Enlaces externos

• Example with VBA code

• Datos: Q1819156